#include <cmath>
#include <thread>
#include <future>
#include <iostream>
#include <vector>
#include <iomanip>
#include <numeric>
#include <cassert>
#include "timer.hpp"

std::mutex io_mutex;

// global
struct sync
{
    sync( std::ostream& os )
        : os(os)
        , lock(io_mutex) {}
    template <class T>
    std::ostream& operator<<(T const& x)
    {
        return os << x;
    }
private:
    std::ostream& os;
    std::lock_guard<std::mutex> lock;
};



/* Riemann computes the riemann sum for a given function f between
 a given 'a' and a given 'b' with a given subinterval length. This
 method is for threaded use, so the result can be accessed from different
 threads. To avoid race conditions, mutex and lock_guard are used.
*/
void riemann(std::pair<long double,
             std::mutex>& result, std::function<double(double)> f,
             double a, double b, double intervalLength)
{
    assert(b>=a);
    assert(intervalLength > 0);

    double xi = a;
    long double sum = 0.;
    sync(std::cout) <<" b is zoveel: "<< b<< " en dx is: "<< intervalLength <<std::endl;
    int nIntervals = static_cast<int>((b-a)/intervalLength);

//This for loop iterates over the intervals. In each step the surface of the interval is added.
    for(int i = 1; i<=nIntervals; ++i)
    {
        sum += f(xi + intervalLength/2)*intervalLength;
        xi +=intervalLength;
    };

    std::lock_guard<std::mutex> l (result.second);
    result.first +=sum;
}



int main()
{
timer t;
t.start();
    double a = 1.;
    double b = 4.;
    double nIntervals = 20.;



// decide how many threads to use
    std::size_t const nthreads = std::max(1u, std::thread::hardware_concurrency());

// Lambda function for the evaluation
    auto f = [=] (double x)
    {
        return log(x)*sqrt(x);
    };

    std::vector<std::thread> threads(nthreads);
// let us just use a single result
    std::pair<long double, std::mutex> result;
    result.first = 0.;
// step is the number of subintervals each thread has to process.
    long double const step = (nIntervals + 0.5l) /nthreads;

    // dx is the length of 1 subinterval.
    double dx = (b-a)/nIntervals;

    for (unsigned i = 0; i < nthreads; ++i)
        threads[i] = std::thread(
                         riemann, std::ref(result),f,
                         a+ i *step*(b-a)/nIntervals,
                         a +  (i+1) * step *(b-a)/nIntervals,dx
                     );

    for (std::thread& t : threads)
        t.join();
t.stop();
// no need to lock here
    std::cout << "The result is " << std::setprecision(18)
              << result.first << " and did it in " << t.get_timing()<< " seconds" <<std::endl;

    return 0;

}



